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ABSTRACT 

We present a fast method for estimating the cosmic microwave background polar- 
ization power spectra using unbiased estimates of heuristically-weighted correlation 



functions. This extends the O(N^) method of Szapudi et al. (2001) to polarized 
■ data. If the sky coverage allows the correlation functions to be estimated over the full 

range of angular separations, they can be inverted directly with integral transforms 
and clean separation of the electric (E) and magnetic (B) modes of polarization is 
obtained exactly in the mean. We assess the level of E-B mixing that arises from 
apodized integral transforms when the correlation function can only be estimated for 
a subset of angular scales, and show that it is significant for small-area observations. 
| We introduce new estimators to deal with this case on the spherical sky that preserve 

, E-B separation; their construction requires an additional integration of the correla- 

tion functions but the computational cost is negligible. We illustrate our methods with 
application to a large-area survey with parameters similar to Planck, and the small- 
area Background Imaging of Cosmic Extragalactic Polarization experiment. In both 
cases we show that the errors on the recovered power spectra are close to theoretical 
expectations. 
^ ■ 

C$ . Key words: cosmic microwave background - methods: analytical: - methods: nu- 
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1 INTRODUCTION stroy polarization on scales that are then sub-Hubble, but 

generates additional large-angle polarization (Zaldarriaga 

With the recent detection of polarization in the cosmic mi- ,„„-> ,. , n . n c , . 

, , . _ , _. , 1997). the potential cosmological returns Irom polariza- 

crowave background CMB by the Degree Angular acale ,. , , . , „ , . ,. , . , T , 

v ' tion observations are nigh. Current polarization data (Ko- 
Interterometer (DASI; Kovac et al. 2002 , and the detec- , „„„„ T , . . , nnno \ , , ,, 

v " vac et al. 2002; Kogut et al. 2003) already allows a stringent 

tion of the temperature-polarization cross correlation by the n , r ,. r „ 

, ; test to be made ol the paradigm lor structure lormation 

Wilkinson Microwave Anisotropy Probe (WMAP) (Kogut r . .,. „ TT , , , . ,. , ,. „ 

, „„„„> , . , r ■ i ■ ■ rrom initially super-Hubble, passive, adiabatic tluctuations. 

et al. 2003), the immediate goal lor upcoming polarization „ ,, , , ,. , ,r rt rA n r i ^ 

' ° ° furthermore, the detection by WMAP ol large-scale power 

experiments is to map accurately the polarization power . , . ,. n 

, f m the temperature-polarization cross power spectrum has 
spectra over a wide range ol angular scales. CMB polariza- . , , . . . .... , 

provided new constraints on the reiomzation process, and 

tion is generated at last scattering irom the local quadrupole , , , . , , „ , , . 

° helped hit major degeneracies that aiiect the determma- 

moment of the photon total intensity. The expected r.m.s. ,. c , . , r ^nirr, . 

J tion ol cosmological parameters Irom the CMB tempera- 
polarization is ~ 6.4 uK, peaking around an angular scale ol . , T n , , , 

, ture anisotropics. (In particular, the degeneracies between 
~ lOarcmin (the angle subtended by the width ol the last .... ,. , , , , , . , . 

v ° J reiomzation optical depth and the scalar spectral index and 

scattering surface), making mapping a challenging prospect. .. .. , ... , % „ , .. , c c 

° ' ° ° i i gravitational wave amplitude.) Potential returns Irom lu- 

Further scattering once the Universe reiomzes tends to de- , . , . . , , ,. s , , ,. c 

° ture, more precise, observations include (l) detection ol the 

clean signature of a stochastic background of gravitational 

waves (Seljak & Zaldarriaga 1997; Kamionkowski, Kosowsky 

* E-mail: gchon@mrao.cam.ac.uk & Stebbins 1997) and hence fine selection between infla- 

1 http://map.gslc.nasa.gov/ 
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tion models (Kinney 1998); and (ii) evidence for weak grav- 
itational lensing through the distortion of CMB polariza- 
tion on small scales (Zaldarriaga & Seljak 1998; Benabed, 
Bernardeau & van Waerbeke 2001). Estimating the polar- 
ization power spectra is an important intermediate step in 
achieving these science goals. This paper addresses that 
problem, presenting a fast, robust method for estimating 
the power spectra from large datasets in the presence of 
real-world complications, such as incomplete sky coverage 
and inhomogeneous noise. 

The accurate analysis of CMB data places strong de- 
mands on the statistical methods employed. Even for to- 
tal intensity data, which is simpler than polarization, the 
extraction of the power spectrum from upcoming mega- 
pixel datasets with standard maximum-likelihood methods 
(e.g. Bond, Jaffe & Knox 1998) is beyond the range of any 
supercomputer. [The operations count scales as the number 
of pixels cubed, JVp ix , while the storage requirements are 
0(iVp ix ).] In the search for fast alternatives to brute-force 
maximum-likelihood power spectrum estimation, two broad 
approaches have emerged. In the first, experiment-specific 
symmetries are exploited to make the brute-force analysis 
tractable, or, if the symmetries are only approximate, to 
pre-condition an iterative solution to the likelihood max- 
imisation. An example of the former is the ingenious "ring- 
torus" method of Wandelt & Hansen (2003), while the lat- 
ter was pioneered by Oh, Spergel & Hinshaw (1999) during 
the development of the pipeline for the WMAP satellite. 
The second class of methods sacrifice optimality in favour 
of speed by adopting a more heuristic weighting of the data 
(such as inverse weighting with the noise variance). An es- 
timate of the underlying power spectrum is then obtained 
from the raw, rotationally-invariant power spectrum of the 
weighted map (the pseudo-Cis; Wandelt, Hivon & Gorski 
2001) either by a direct linear inversion (Szapudi, Prunet & 
Colombi 2001; Hivon et al. 2002) or with likelihood methods 
(Wandelt et al. 2001; Hansen, Gorski & Hivon 2002). The 
linear inversion, which yields estimators quadratic in the 
data, can be performed directly in harmonic space (Hivon 
et al. 2002), or, more simply, by first transforming to real 
space (i.e. by constructing the correlation function) and then 
recovering the power spectrum with a (suitably apodized) 
integral transform (Szapudi et al. 2001). The estimation 
of polarization power spectra is less well explored than for 
total intensity, although all of the above methods can, in 
principle, be extended to handle polarization. To date, only 
brute-force maximum likelihood (Kovac et al. 2002, Munshi 
et al. in preparation), minimum- variance quadratic estima- 
tors (Tegmark & de Oliveira-Costa 2002), pseudo-Cj meth- 
ods with statistical (Hansen & Gorski 2003) or direct inver- 
sion (Kogut et al. 2003) in harmonic space, and real-space 
correlation function methods (Sbarra et al. 2003) have been 
demonstrated on polarized data. Of these, only the pseudo- 
C'i methods are fast enough to apply many times (e.g. in 
Monte-Carlo simulations) to mega-pixel maps. In this paper 
we extend the fast correlation-function approach of Szapudi 
et al. (2001) to polarized data. 

A new problem that arises when analysing polarized 
data, that is absent for total intensity, is the decomposi- 
tion of the polarized field into its electric (E; sometimes 
denoted gradient) and magnetic (B; alternatively curl) com- 
ponents. The scientific importance of this decomposition is 



that primordial magnetic polarization is not generated by 
density perturbations, and so in standard models is sourced 
only by gravitational waves (Seljak & Zaldarriaga 1997; 
Kamionkowski et al. 1997). However, with incomplete sky 
coverage, separating the polarization field into electric and 
magnetic components is no longer straightforward. Exquisite 
monitoring of leakage between E and B in analysis pipelines 
will be required if primordial B polarization is to be de- 
tected down to the fundamental confusion limit set by cos- 
mic shear (Kesden, Cooray & Kamionkowski 2002). The 
question of performing the E-B separation on an incomplete 
sky has received considerable attention recently (Zaldarriaga 
2001; Lewis, Challinor & Turok 2002; Bunn 2002; Chiueh & 
Ma 2002; Bunn ct al. 2003). Methods are now available for 
extracting pure measures of the E and B fields which can 
then be used for subsequent power spectrum estimation. An 
alternative approach is to perform a joint (i.e. E and B) 
power spectrum analysis of the original polarization data, 
removing the need for an additional stage in the analysis 
pipeline and the non-optimalities that this may introduce. 
The efficacy of maximum-likelihood methods for performing 
the E-B separation is explored by Munshi et al. (in prepa- 
ration). However, the computational demands of likelihood 
methods, and the difficulty in monitoring E-B leakage in a 
non-linear analysis, motivates the development of fast, unbi- 
ased methods. The correlation-function based approach we 
develop here, motivated by Crittenden et al. (2002), has a 
significant feature in that, with a little post-processing of 
the correlation functions, leakage between E and B can be 
eliminated in the mean, even for observations covering only 
a small part of the sky. The separation is exact in the mean. 

The outline of this paper is as follows. In Sec. 2 we 
review the polarization functions on the sphere and their 
relation with the power spectra. Section 3 presents a fast, 
0(Np{^) method for computing unbiased estimates of the 
correlation functions allowing for heuristic weighting of the 
data, and describes power spectrum recovery for large-area 
surveys where the correlation function can be estimated for 
the full range of angular separations. We illustrate our meth- 
ods by applying them to a survey mission with similar pa- 
rameters to those for Planck 2 . In Sec. 4 we provide a careful 
analysis of the effect of incomplete coverage of the correla- 
tion functions on the direct extraction of the power spectra 
with apodized integral transforms. By constructing the rel- 
evant window functions for small-area observations we show 
that leakage from E to B can be a significant problem. We 
remedy this deficiency of the method in Sec. 5, where we 
construct functions from integrals of the original correla- 
tion functions that contain signal contributions from only 
E or B in the mean. These functions can be safely inverted 
with apodized integral transforms to obtain properly sepa- 
rated estimates of the E and B power spectra. We apply 
this new estimator to a model of the Background Imaging 
of Cosmic Extragalactic Polarization experiment (BICEP) 3 
experiment, and show that it produces error bars close to 
the theoretical expectations. Our conclusions are given in 
Sec. 6, and the Appendix contains some technical results 
on the analytic normalisation of the correlation functions 



2 http:/ /sci.esa.int/homc/planck/ 

3 http://bicep.caltech.edu 



© 2003 RAS, MNRAS 000, 1-14 



Fast polarization power spectrum estimation 3 



estimators for uniform weighting on azimuthally-symmetric 
patches. 

Throughout this paper we illustrate our results with a 
flat, ACDM model with concordance parameters Q^li 2 = 
0.022, Q c h 2 = 0.12, fi A = 0.7 giving the Hubble constant 
h — 0.69. The primordial scalar curvature and tensor spectra 
are scale-invariant and have ratio r = 0.31 (making the ratio 
of tensor to scalar power in temperature anisotropies no = 
0.16 at I = 10) 4 . We ignore the effects of weak gravitational 
lcnsing. We consider two models for the ionization history: 
no reionization and full reionization at redshift z IC — 6. Note 
that had we adopted a model with earlier reionization, e.g. 
z rc ~ 15 as favoured by WMAP data (Kogut et al. 2003), 
the problem of E-B mixing described in Section 4 would 
have been further exacerbated by the additional large-scale 
E power. 



2 POLARIZATION CORRELATION 
FUNCTIONS ON THE SPHERE 

Stokes parameters Q{h) and U{h) are defined for a line of 
sight n with the local ir-axis generated by 6 and the local y- 
axis by —<j>. Here 9 and 4> are the basis vectors of a spherical- 
polar coordinate system. A right-handed basis is completed 
by the addition of the radiation propagation direction — n. 
The polarization P = Q + iU is spin —2 (Newman & Penrose 
1966) and can be expanded in spin-±2 harmonics as (Seljak 
& Zaldarriaga 1997) 



(Q ± iU)(n) = ^{Ei m T iBi m ) T 2Yi m {h). 



(1) 



Reality of Q and U demands E\ m = (— l) m i5;_ m with 
an equivalent result for Bi m . Under parity transformations, 
(Q±iU)(h) -> {Q + iU){-h) so that E im has parity (-1)' 
(electric), but Bi m has parity (— 1)' +1 (magnetic). The tem- 
perature is a scalar field and so can be expanded in spherical 
harmonics with multipoles T; m . In an isotropic- and parity- 
invariant ensemble the non-vanishing elements of the polar- 
ization covariance structure are 



(Ei m E Vrn i) 
(Bi m B*, m ,) 



(2) 
(3) 
(4) 



If the direction hi corresponds to angular coordinates 
{0i, cj>i), and similarly for n 2 , then the SO(3) composition 



D- 1 {<t>i,0i,O)D{<f> 2 ,02, 0) = D(a, (3, - 7 ) 



(5) 



determines (3 (0 < /3 < ir), the angle between hi and A2, a, 
the angle required to rotate 0{hi) in a right-handed sense 
about hi onto the tangent (at hi) to the geodesic connecting 
hi and n 2 , and 7, defined in the same manner as a but 
at tj.2- Making use of the relation between the Wigner-D 
matrices (e.g. Varshalovich et al. 1988) and the spin-weight 
spherical harmonics, 



6,0,0) = (-l)" 



47T 

21 + 1 



sYi m {h), 



(6) 



4 We define r as the ratio of the amplitudes As and At of the 
primordial curvature and tensor power spectra, following the con- 
ventions of Martin & Schwarz (2000). 



we obtain the following representation of equation (5): 
D l ss ,{a,P,—y) = ^ 21 + 1 = Y i*m{ni)s'Yim{h 2 ). 



(7) 



With this result, the two-point correlation functions for lin- 
ear polarization evaluate to (Ng & Liu 1999) 

(P{hi)P{h 2 )) = ^^±l (C f-Cf)4- 2 (/3), (8) 
(P*{hi)P{h 2 )) = ]T^±i(cf +Cf)4 2 (/3), (9) 



(T(m)p(n 2 )) = j2^cr E d i 20 m 



(10) 



where d„ 



are the reduced D-matrices. Note that 



(T(ni)P(n 2 )> = (P{hi)T{h 2 )). The quantities 

Pirn) = e 2ia P{h!), (11) 

P(n 2 ) = e 2i ^P(n 2 ), (12) 



are the polarizations defined on local bases with the in- 
direction along the geodesic between hi and n 2 . With these 
rotations, the correlation functions depend only on the angle 
[3 between the two points. Note that (P*(ni)P(n 2 )) is real, 
which follows from statistical isotropy, while (P(ni)P(n 2 )) 
and (T(ni)P(n 2 )) are only real if the universe is parity- 
invariant in the mean. In the presence of parity violations, 

2/ + 1 . 



<P(m)P(n 2 )> 



(T(m)P(n 2 )) = 
where 

{El m B* m ) — Su'Sr, 



1 



47T 



xd2-a(/?)] 
21 + 1 



(13) 



E 



r ,t 



47T 



iCT B )d l 20 {!3), (14) 



(TlmBi m ) = du'Smm'Cl .(15) 

The correlation functions of the (rotated) Stokes parameters 
can be found directly from those for P. Defining 



= (p{hi)p{h 2 )), 

= (P*{hi)P{h 2 )), 

ix{(3) = (T(m)P(n 2 )), 



we have 

(Q{hi)Q{h 2 )) = 

(U{hi)U{h 2 )) = 

(Q{hi)U{h 2 )) = 

(T{hi)Q{h 2 )) = 

{T{hi)U{h 2 )) = 



i[C+(/3)-^- (/?)], 



(16) 
(17) 
(18) 

(19) 
(20) 

(21) 

(22) 
(23) 



In Fig. 1 we plot the correlation functions £±(/3) and £x{/3) 
for the cosmological models described in Sec. 1. The damp- 
ing of the polarization power on linear scales that are sub- 
Hubble at the epoch of reionization (in this case at z = 6) 
is just discernible in the correlation functions in Fig. 1. The 
additional large-scale power due to reionization makes a neg- 
ligible contribution to the correlation functions for the an- 
gular range we have plotted. 
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z r =0.0 
z =6.0 




Figure 1. The correlation functions £+(/?) (top panel), £-((3) 
(middle), and £x(P) (bottom). The cosmological model is as de- 
scribed in Sec. 1; the solid lines are for no reionization, while the 
dashed have complete reionization with z Ie = 6. The angle fi is 
in degrees. 



We can invert equations (9), (13) and (14) using the 



orthogonality of the d„ 



dLm' W)d l mm > (P) d COS [3 = ^ — - S W , 



(24) 



to determine the power spectra from the correlation func- 
tions: 

Cf - Cf - 2iCf B = 2tt y ((3)d l 2 _ 2 (f3) d cos /?, (25) 

i 



Cf + C? = 2tt / C+(/3)4 2 (/3)dcos/3, (26) 



C, + iCi 



2n / Cx(/3)4o(/3)dcos/3. (27) 



Since those reduced D-matrices that appear in these equa- 
tions are are polynomials in cos /3 for, the integrals can 
be performed essentially exactly for band-limited data by 
Gauss-Legendre quadrature. 



3 FAST CORRELATION FUNCTION 
ESTIMATORS 

If we had available unbiased estimates of the various cor- 
relation functions for all angles ^ f3 ^ n, we could ob- 
tain unbiased estimates of the power spectra by performing 
the inversions in equations (25)-(27). Direct evaluation of 
the correlation functions (e.g. Sbarra et al. 2003) requires 
0(iVp ix ) evaluations and is complicated by the need to per- 
form a rotation to the appropriate basis for each pair of 
points. Here we consider an O(N^) method based on fast 
spherical transforms. This method generalises that of Sza- 
pudi et al. (2001) to polarization fields. 

We consider an arbitrary weighting of the noisy polar- 
ization field P(h) and the noisy temperature field T(h) with 
some weight functions wp(n) and wr(fi) respectively. The 
weight is zero for those pixels in regions that are either not 
observed or are removed from the map due to foreground 



contamination. In this paper we only consider real weight- 
ing of the polarization field, thus preserving the direction of 
polarization at any point; relaxing this condition is straight- 
forward if required. In the presence of instrument noise, the 
weights allow for a heuristic pixel-noise weighting of the 
data. We start with the following estimators for the signal- 
plus-noise correlations: 



C+(ijj) — Ap(ip) J dnidn 2 [S(ni ■ n 2 — cosi/>) 

x wp(ni)wp(n 2 )P*(ni)P(n 2 )l{28) 

C-(ip) — Ap{ip) J dnidn 2 [8{h\ ■ n 2 — cosi/>) 

x w P (ni)wp(n 2 )P(rii)P(n 2 )], (29) 

Cx(4>) — Ax(4>) J dnidn 2 [5(ni ■ n 2 — cosi/>) 

x w T (ni)w P (n 2 )T{ni)P{n 2 )]. (30) 

The delta functions ensure that we only consider those 
points that have angular separation tp. The normalisations 
A(ip) are chosen so that our correlation function estimators 
are unbiased in the absence of noise. This requires 



A P (ij) 



dnidn 2 [5(ni ■ n 2 — cosV>) 
x wp(ni)wp(n 2 )], 

dnidn 2 [5(n-i ■ n 2 — cosV>) 

X WT(ni)wp(fl2)]- 



(31) 



(32) 



These expressions for the correlation functions and normal- 
isation factor can be simplified by using the completeness 
relation 

2/ -(- 1 

^ 2 d l mn (f3)d l mn {ip) = S(cosf3 - costp) (33) 

max( |m | , | n| ) 

to substitute for the delta functions. To evaluate C+(ip) we 
set m — n — 2, so the integrand in equation (28) involves 

wp(ni)P* ' (hi)d l 22 (/3)P(n 2 )wp(h 2 ) 

= P*(n 1 )D l 22 (a,f3,-j)P(n 2 ), (34) 

where cos/3 = fi\ • n 2 , and we have used equations (11) 
and (12). Here, P(n) = wp(n)P(n) is the weighted polar- 
ization field on the (polar-)coordinate basis. We can now 
use equation (7) to express the D-matrix in terms of spin- 
weight harmonics. Performing the angular integrals extracts 
the spin-weight 2 (pseudo-)multipoles of the weighted, noisy 
polarization field, defined by 



P(n) = y^(g; m - iB lm )- 2 Y lm (n), 

lm 

P*(n) = J2(E lm + iB lm ) +2 Y lm (n), 



leaving 

C+{i>) = 2nA P (i>) ^ d 22 {iP)\Ei m + iB lm \ 2 . 



Introducing the real pseudo-C;s for the weighted fields: 



Cf = 



^ \Ei m \' 2 , 



2l + l^' 



(35) 
(36) 

(37) 

(38) 
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c 



c, 



m 

^— j" Ei m Bt m = ^ - - ^ BlmEtm, (40) 

m m 



21 + 
21 



21 

we can write 



(41) 



(42) 



To evaluate C- (tp) we follow the same procedure, but with 
m = —2 and n = 2 in equation (33). The result is 

C-(^) = 27rAp(^)^(2/+l)4- 2 (V')(Cf-Q B -2iCf B ).(43) 
i 

Finally, for Cx (tp) we take m = 2 and n = to find 

Cx(i>) = 2nA x (<p)J2(2l + Mo(4>)(C? E -iCj B ). (44) 
i 

The normalisation factors A(i/>) can be evaluated by taking 
m = n — in equation (33), e.g. 



Ap(V) 

where 



■■ 2tt^(2Z + l)Pi(cosV)wp,i, 



2/ 



— y 

+i ^ 



(45) 



(46) 



with u>p,; m the (spin-0) spherical multipoles of the weight 
function wp(n). Note that we have used d l 00 (tp) — P;(cosi/>) 
where Pi (x) is a Legendre polynomial. Once the mean noise 
contribution (noise bias) is removed from the estimators 
C(ip) [leaving unbiased estimators of the signal correlation 
functions £(ip)], we can use equations (25)-(27) to compute 
estimates of the power spectra. The real parts of £(?/>) give 
estimates of Cf , Cf and C; X , while the imaginary parts of 
£-(i/>) and £,x(ip) can be used to estimate Cf B and Cj B 
and hence test for parity violations. 

The full set of pseudo-multipoles can be obtained effi- 
ciently in 0(Np{^ log iVpi x ) operations using fast spherical 
transforms such as those implemented in the HEALPix 5 and 
IGLOO(Crittenden & Turok 1998) packages. (Our current 
implementation employs HEALPix.) To remove the noise bias 
from C(ip) it is generally most efficient to resort to Monte- 
Carlo simulations of pure noise fields (Szapudi et al. 2001). 
(An exception is the case where the noise is uncorrelated 
between pixels; see below for details.) The ensemble mean 
of these pure-noise correlation functions can be subtracted 
from C(t/j) to yield (asymptotically) unbiased estimates of 
the signal correlation functions. Monte-Carlo estimation of 
the noise bias provides a robust means of dealing with dis- 
cretisation effects due to the chosen pixelisation. Monte- 
Carlo methods also offer the simplest method of computing 
the variance of the power spectrum estimates. In the pres- 
ence of uncorrelated noise it is straightforward to proceed 
analytically with the noise contribution to the variance, but 



http: / /www. eso.org/science/healpix/ 



the cosmic variance contribution is complicated by the pres- 
ence of signal correlations. 

For the simple case of noise that is uncorrelated between 
pixels it is straightforward to compute the noise bias ana- 
lytically. For simplicity consider noise that is uncorrelated 
between Q, U and T, and has equal variance in Q and U. If 
the noise variance of the Stokes parameters per solid angle is 
ap(n p ), then in the continuum limit the polarization noise 
correlations can be summarised by 



<iMni)fft(n„)) 
(Pjv(ni)Pjv(n 2 )) 
(T^(ni)Pv(n 2 )) 



2a P S(ni — n 2 ) 

0, 

0, 



(47) 
(48) 
(49) 



where Pjv(n) is the spin —2 noise, and Tjv(n) is the noise on 
the temperature. As the noise is uncorrelated between pixels 
its mean effect on correlation function estimates is confined 
to zero separation: 



<AC+(V0> 



<ACx(V)> 



A P (0)5(l-cosV) 



• / (ln%(n)2ffp(n), 



0, 
0, 



with 
1 

aRo) 



= 2tt dn w P (n). 



(50) 

(51) 
(52) 



(53) 



Here, (AC) is the mean noise contribution to the estimators 
C. Making use of d' mm , (0) = S mm i we find that the non- 
zero noise biases in the estimates of the power spectra in 
the continuum limit are 



{ ACf) = (ACf)= Idflw2p{fl)aUfl) 



J dn Wp(n) 



(54) 



However, we would recommend removing the noise bias with 
Monte-Carlo techniques even for simple, uncorrelated noise. 
This is to ensure that the effective band limit introduced on 
the noise by computing the correlation functions via pseudo- 
Cis up to some finite Z max is properly accounted for. 

3.1 Application to large-area surveys 

As an application of our method we consider extracting the 
power spectra from simulated maps obtained with a full-sky 
survey with pixel noise similar to that expected for Planck. 
To be specific, we assumed uncorrelated pixel noise on Q 
and U with r.m.s. 6.95 fiK in a 10-arcmin by 10-arcmin 
pixel. We adopted a beam size of 10 arcmin, somewhat 
larger than the polarization-sensitive channels of the Planck 
High- Frequency Instrument, so to ensure oversampling of 
the beam at HEALPix resolution iVside = 1024. We ignored 
the variation in pixel noise across the map, but this could 
easily be included in our simulations at no additional com- 
putational cost. Noise correlations could also be included 
easily if fast simulation of noise realisations were possible. 
We made a constant-latitude Galactic cut of ±20°. The un- 
derlying cosmological model was as described in Sec. 1 and 
we assumed reionization at z = 6. We adopted a uniform 
weighting scheme motivated by the constant variance of the 
noise. 
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Figure 2. Recovered Cf power spectrum for a large-area survey 
with a ±20° Galactic cut. The points are flat band-power esti- 
mates, with Al = 10, from a single simulation; the shaded region 
shows the ±<r region on the basis of equation (55). 

The recovered power spectrum Cf is shown in Fig. 2. 
We computed estimates of the correlation function at the 
roots of a Legendre polynomial from the pseudo-C;s (ob- 
tained with the fast spherical transforms in HEALPix). Inver- 
sion of the correlation functions was performed with Gauss- 
Legendre integration. We averaged the recovered 1(1 + 1)C; 
to form flat band-power estimates with a Al = 10, and the 
results of one simulation are shown as the points in Fig. 2. 
We adopted a Al = 10 to ensure the errors were essentially 
uncorrelated. The shaded area in Fig. 2, centred on the true 
spectrum smoothed with a 10-arcmin beam, encloses the 
±<r error region based on the rule of thumb (generalised 
from Hivon et al. 2002 for the temperature case) 

ACf « ^(Cf +Ni). (55) 

Here Ni is the (full-sky) noise power spectrum, and vi = 
Al(2l + l)/ s kyW2/u)4 is the effective number of degrees of 
freedom in a band of width Al on a fraction of the sky 
/sky, where 47r«;;/ s ky = J Wp(n) dn. Note that vi takes no 
account of the loss of degrees of freedom associated with 
disentangling E and B polarization, since the fractional loss 
of modes is small for / s k y close to unity (Lewis et al. 2002). 
The scatter of points in the simulation is broadly consistent 
with that expected on the basis of the theoretical errors. A 
more detailed analysis of the optimality of our method must 
await comparison with optimal, maximum-likelihood codes 
when these become available at sufficiently high resolution. 



a 90° band about the Galactic plane. In the case where the 
correlation functions cannot be estimated for all separations 
ip, estimating the power spectra by direct integration (e.g. 
equations 25 and 26) over the observed range will introduce 
window functions ± 2 Kui such that 

(Cf ± Cf) = J2 ±* K u' ( C v ±Cf). (56) 
v 

Fourier ringing can be reduced by pre-multiplying the corre- 
lation functions with a scalar apodizing function f(ip) prior 
to integration, in which case the window functions take the 
form 

±2 K W = —Z- j /(/3)4± 2 (/3)4± 2 (/3)dcos/3, (57) 

where the integral is over the range of angles for which the 
correlation functions can be estimated. Note that the win- 
dow functions are not symmetric but rather satisfy 

{2l + l) ±2 K w = (2l' + l) ±2 K in . (58) 

Introducing the sum and difference window functions, 
±K H / = (2K111 ± - 2 K H r)/2, the means of the estimated 
power spectra are related to the true spectra by 

(Cf) = Y,{+K tt ,Ci? + -K u ,Cf), (59) 
v 

(Cf) = ^(-K w Cf + +K w Cf). (60) 
v 

The window function -K H i controls the mixing of E and B 
polarization. Recent results from DASI (Kovac et al. 2002) 
are in line with theoretical expectations that B polariza- 
tion should be sub-dominant, so that cross contamination 
due to partial sky effects is proportionately more troubling 
for B polarization than for E. While not presenting a fun- 
damental problem for cosmological parameter extraction, a 
non-zero -Kui makes interpretation (and presentation) of 
the estimated Cf awkward. Mixing can obviously be elim- 
inated by pre-multiplying the estimates Cf ± Cf with the 
inverse of the window functions ± 2 K a / , but this inversion is 
awkward in practice due to the ill-conditioned nature of the 
window functions when coverage of the correlation functions 
is incomplete. In Section 5 we introduce a simple, robust 
technique for extracting the power spectra from correlation 
functions which eliminates mixing in the mean (i.e. produces 
a zero -K u i). Before turning to that, in the following subsec- 
tions we first explore the properties of the window functions 
given by equation (57) , and the circumstances under which 
mixing is significant. 



4 WINDOW FUNCTIONS FROM 

INCOMPLETE CORRELATION FUNCTIONS 

In this section we construct the window functions that arise 
when the correlation functions can only be estimated over a 
limited angular range. We shall concentrate on the polariza- 
tion auto-correlations; the generalisation to the polarization- 
temperature cross-correlation is straightforward. 

If unbiased estimates of the correlation functions are 
available over the full angular range (0, 7r), they can easily be 
inverted to obtain unbiased power spectra. This case would 
describe full-sky experiments with a cut excising less than 



4.1 General properties 

If we define the apodizing function f(tp) to be zero outside 
the observed range of if), and perform a Legendre expansion 

the window function 2 K u i reduces to 

^^E^+^O - 2 o) 2 ' (62) 
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while -iK u t has an additional factor of (— +V) in the 
summation. The array in brackets in equation (62) is a 
Wigner 3-j symbol arising from the integral of a product 
of three rotation matrices. If the apodizing function is ef- 
fectively band limited to L max , the window functions vanish 
for |2 - l'\ > i max . 

The normalisation y] f , ±2-K«' of the window functions 
is also of some interest. For 2K H i we can perform the sum 
over I' in equation (62) directly using the orthogonality of 
the 3-j symbols (Varshalovich et al. 1988) to find 



E 2 ^' = E 



2L + 1 



h = /(o). 



(63) 



The last equality follows from Pi(l) = 1, or, more directly, 
by employing + l/2)d 2 2 (P) = <5(cos/3 — 1) in equa- 

tion (57). The normalisation of -2-K)i' is a little more in- 
volved. We start with the result 



E ^7^4 - 2 (/3) = <5(cos/3 - 1) + csc 2 (/3/2), 



(64) 



which follows by summing equation (89) of Section 5 over I' . 
If we now sum equation (57) over I' , and use equation (64), 
we find that 



V J 



f((3) csc 2 (/3/2)d 2 - 2 (/3)dcos/3, 



(65) 



where we have used ^2-2(0) = [and the assumed regular- 
ity of f(P)]- The function csc 2 (/3/2)d 2 -2(P) ls a polynomial 
in cos p and so the integral can easily be evaluated numeri- 
cally by e.g. Gauss-Legendre integration for smooth apodiz- 
ing functions. To make further progress analytically we in- 
sert the Legendre expansion of /(/?) in equation (65) and 
use the differential representation of the reduced D-matrices 
(e.g. Section 4.3.2 of Varshalovich et al. 1988). Repeated in- 
tegration by parts then establishes the result 



E 



- E 



2L + 1 



h 1 - 4 



L(L + 1) 
1(1 + 1) 



+3 



(L + 2)! (Z-2)! 
{L-2)\ (1 + 2)} 



(66) 



If f(P) is effectively band-limited, for I 3> L max we have 
y^ ; , -2K111 ~ /(0). In this limit, the normalisation of -Kw 
is much smaller than that of +K H /, and mixing of E and 
B-power is suppressed in the mean (Bunn 2002). 



4.2 No apodization 

Consider the case where the correlation functions can be 
estimated in the range (0, Anax). If we apply no apodization 
to the correlation functions, we obtain window functions 



21' 



±2VV U 



4s ±2 (0)4 ±2 (/?)d cos/ 



(67) 



For I / I' this integral can be evaluated directly since the 
d l mn are eigenfunctions of a self-adjoint operator. The result 



Ww = 



21' + 1 



±2Wll 



COS 2 Pn 



2 1(1 + 1) - I' (I' + 1) \ dcos/3 



dd 2 ±2 



dd 2 ±2 

" J S"2±2 

dcos p 



(68) 



For / = V the integral can be evaluated recursively as de- 
scribed in Appendix C of Lewis et al. (2002). The window 
function - W u i = (2W111 — -2W H i)/2 can be evaluated di- 
rectly for all I and I' (Lewis et al. 2002): 



-Wn 



21' + 1 



(UIU V + VlVl>)\f 



where the vectors 



ui(f3) = 
vi(/3) = 



0-2)! . d 
sin p - 



(I + 2)\ 



d/3 \sm(3 



(I ~ 2)! V3 j 
(i + 2)!sin/3 2( 



(P)- 



(69) 

(70) 

(71) 
= 



Both vectors vanish for /3 max = 7r to ensure that - W H 
when the full angular range (0, 7r) is considered. 

Some representative rows of the window functions ±Ww 
are shown in Fig. 3 for /3 m ax = 20° (corresponding to e.g. 
observations over a circular patch of radius 10°). Note that 
+ Whi is localised around I = I' (with width varying inversely 
with /3 max ), while — Ww shows no localisation. Equation (69) 
shows that, considered as a matrix, - W H i is of rank 2, so 
the rows of the window function are constructed from lin- 
ear combinations of up and vy . The approximate scaling of 
-Ww with I for fixed I' , that is evident in Fig. 3, arises be- 
cause the vector u t oscillates with larger amplitude than vi 

for / > 1/Anax- 

In Fig. 3 we also show the mean of the estimated power 
spectrum (Cj 3 ) obtained by multiplying the window func- 



cosmological models detailed in Sec. 1. The true spectra are 
convolved with a Gaussian beam of full- width 10 arcmin 
at half-maximum. In the case of no reionization, the mean 
of the recovered Cf is a faithful representation of the true 
spectra. This is because (i) the inner products between Cfi 
and either of (/' + l/2)u;/ and (I' + l/2)v t i are sufficiently 
small that the leakage from E polarization causes only a 
small amplitude oscillation in the recovered Cf ; and (ii) for 
+ Wui sufficiently localised compared to the scale of features 
in Cy , we can approximate their product by 



E 



E 



Wn 



Cf (/>1/Anax). (72) 



For the second approximation note that the window function 
+ Wni inherits its normalisation from that of 2W111 (which is 
unity) and -2W1V . For the latter we use equation (65) and 
the differential representation of d 2 _ 2 to find 



E- aW « 

V 

X \^Jl(l + l)d 



1 + 2 



0-2)! 
(2 + 2)! 



COt(/3 ma x/2) 



\_ 2 (P)-cot(p/2)d l _ 2 (P) 



(73) 



For large I ^> l//3 max we find Y]., -2Ww ~ 1 (see also the 
discussion after equation 66). Reionized models are more 
problematic since they have additional large-scale power in 
E polarization (and so are more sensitive to the truncation 
of the correlation functions at /3 m ax). The effect of this large- 
scale power can clearly be seen in Fig. 3 for the model with 
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Figure 3. Left: representative rows of the window functions + (top panel) and —Wui (bottom panel) when the correlation functions 
are known in the angular range (0,20°), and no apodization is applied. The solid lines are for I = 30, the dashed for / = 120 and the 
dotted for I = 210. Right: mean recovered Cf (solid lines), obtained from the convolution in equation (60) with the windows ±W;;/ , 
compared to the true Cf (dashed lines). The top panel has no reionization while the bottom panel is a model with full reionization at 
z = 6. 



reionization at z — 6. The large amplitude oscillations in the 
recovered Cf trace those of the vector U!(/3 max ) at large I. 

4.3 Gaussian apodization 

The Fourier ringing evident in the window functions in Fig. 3 
can be reduced by apodizing the correlation functions. Here 
we consider Gaussian apodizing functions, i.e. 

/(/?)=e -/3W). (74) 

The half-width at half-maximum is aV 2 In 2 which should 
be small compared to the cut off /3 max in the correlation 
functions for effective apodizing. The window functions, ac- 
counting for apodization and the finite range (0, /3 ma x) of 
the observed correlation functions, can be written as matrix 
products: 

±iK w = ^2±2F lL±2 W L v, (75) 

L 

where 

11' 4- 1 I' 1 i 
±a F u , = —^—j /(/3)d< ±2 03)4 ±2 03) d cos /3. (76) 

Note that the full window functions ±2K H > are insensitive 
to the behaviour of the apodizing function for j3 > /3 max . 
Note also that the order of the matrix product in equa- 
tion (75) is irrelevant since the window functions commute. 
If the apodizing function is narrow compared to /3 max we 
expect ±iK lv w ±iF lv . 

In Fig. 4 we show representative rows of the sum and dif- 
ference window functions, ±K u i , for /3 max = 20° and a Gaus- 
sian apodizing function with half-width at half-maximum 
equal to /3 max /2 [so /(/3 max ) = 1/16]. These are well ap- 
proximated by Gaussians centred on I — I' with width 1/er. 
The amplitude of the difference window functions are much 
smaller than those of the sum, and for large I the ratio of 
amplitudes oc I /(la) 2 . 

To understand this behaviour, consider the limit where 
(r«l and la ^> 1. In this flat-sky limit we can approximate 



the reduced D-matrices by Bessel functions (Varshalovich et 
al. 1988): 

4 2 (/3) w Jo (IP), 4- 2 (/3)~ J 4 (//3), (77) 
and the window functions ±2-Fh' are easily computed to be 

2 F W » l'a 2 e- a2 ^ +l ' 2) ' 2 h{ll'a 2 ), (78) 

- 2 F U , « l'a 2 e- a2{l2+l ' 2)/2 h(ll'a 2 ). (79) 

The leading-order asymptotic expansions of ±F u i (11' a 2 
1) then follow: 

which reproduce the behaviour seen in Fig. 4. The win- 
dow function 2 Fiy is normalised to unity by virtue of equa- 
tion (63). The normalisation of -iFui can be calculated for 
all I in terms of modified spherical Bessel functions by ap- 
proximating f([3) w exp[— (1 — cosf3)/a 2 ] in equation (65). 
However, the result is cumbersome so we shall only give its 
asymptotic form here (valid for a <ti 1 and la ^> 1): 

+ -.-*■'■). m 

(This result also follows directly from integrating equa- 
tion 79 over I'.) Asymptotically, we then have y] ?< -F H i ~ 
4/ (la) 2 , consistent with Fig. 3. 

We also show in Fig. 4 the mean of the recovered Cf 
for a range of I values with Gaussian apodizing (corre- 
sponding to the window functions in the left-hand panels 
of Fig. 4). The spacing of points is chosen to reflect the l- 
range over which recovered power spectra should be roughly 
decorrelated. While apodizing has clearly removed the high- 
frequency oscillations in the mean of the recovered Cf , it 
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Figure 4. Left: representative rows of ±Kui for /3 max = 20° and Gaussian apodization with half-width at half-maximum equal to 
Anax/2. Right: mean recovered Cf (points) compared to the true Cj 3 (lines) with (bottom) and without (top) reionization. 



has done so at the expense of introducing considerable bias 
due to leakage from E polarization. As expected, apodiza- 
tion has reduced the sensitivity of the recovered Cf to the 
level of large-scale power in E polarization (and hence reion- 
ization; cf. Fig. 3), but has replaced it with a local bias 
w 4Cf/(la) 2 . The bias becomes non-local in models with 
sufficiently early reionization, where the level of large-scale 
E power can be such that it is transmitted to the recov- 
ered B power spectrum for all I through the low-? tail of the 
window function. 

A more effective way to reduce oscillations in Cf with- 
out introducing additional bias due to E-B leakage is to 
recover the power spectra with no apodization, and then 
to post-convolve with a suitably-wide smoothing function. 
Such an approach is illustrated in Fig. 5, where we have 
post-convolved the results in Fig. 3 with the asymptotic 
form of +F u i (equation 80). The Gaussian smoothing pro- 
duces well-localised +K u i window functions, removing the 
Fourier ringing from +Wu>, and reduces the amplitude of 
-Wu> by a factor ~ 40. Only a low- frequency, oscillatory 
bias remains in the recovered Cf in the model with reion- 
ization. Whether this bias is significant depends on the level 
of reionization (and primordial B polarization). 



5 REMOVING E-B LEAKAGE 

Although we are able to reduce the level of cross- 
contamination in the recovered power spectra by simply 
post-convolving them with a suitably-wide smoothing func- 
tion, it is actually straightforward to remove this E-B mix- 
ing exactly in the mean. Crittenden et al. (2002) showed how 
to construct correlation functions on small patches of the sky 
that contain only E or only B modes in the mean. (Their 
work was in the context of weak gravitational lensing, but 
their results are equally applicable to CMB polarization.) In 
this section we extend the central result of Crittenden et al. 
to the sphere, so we are able to handle large-angle polariza- 
tion signals, and demonstrate the methods with simulations 
for an experiment similar to BICEP. 

We begin by considering the function 



(83) 



If we had access to over some range of scales we could 
combine with the real part of £-(/3) to extract the function 



(84) 



which depends only on B-polarization. We could thus re- 
cover an estimate of Cf by integrating unbiased estimates 
°f ~ £- (P) against d l 2 _ 2 [with appropriate apodization 
/(/?)]: 

Cf = ^ J - 3^-(/3)]/(/3)4- 2 (/3)dcos/3. (85) 

Such an estimate would contain no contamination from E 
polarization in the mean (i.e. the difference window function 
-Kui would vanish for arbitrary apodization). A similarly 
unbiased estimate of Cf could be obtained by considering 

£(/?)+»£-(/?). 

The result we now prove is that the function £ (/?) can be 
obtained in the range (0, /3 ma x) from the correlation function 
£+(/?) in the same range by quadrature. We start by inserting 
equation (26) into the summand of equation (83) which gives 



J-i , 



21 + 1 , 
— o — "2 - 



2 (/?)<& (/?')■ (86) 



Our strategy for simplifying the summation in this equation 
is to express d l 2 -2{P) m terms of integrals involving d l 2 2(P), 
and then perform the summation with the completeness re- 
lation, equation (33). Making repeated use of the recursion 
relation (Varshalovich et al. 1988) 

— m + m' cos/3 t 



sin/3 -d l mm i(/3) = ±y/(l + m')(l -m' + 1)4, 
+ 1 y/(l-m>)(l + m' + l)<L m , +1 (/?) 



(87) 



and the relation 

fR s m - m! cc 



= -^{l-m'){l + m' + l)d\ 



mm' + 1 



(88) 
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Figure 5. Window functions (left) and mean of recovered Cf (right) obtained by post-convolving the results in Fig. 3 with the Gaussian 
asymptotic approximation to +Fj;/. 



we find that 



c/3 



A-m = d l 22 {P) - 2(2 y° S 9 f / tan 3 (/3'/2)4 2 (/3W 
sm (p/z) J 



+ 



sec 3 (/3'/2) sin(/?72)4 2 (/3') d/3'. (89) 



Multiplying with (I + l/2)d 2 2(/?') and summing over I we 
find 

m = £+(/?)+ ~ 27^ / dcos/3'C+(/3')sec 4 (/372) 



sm 2 (/3/2) 

2(2 + cos /3) Z* 1 
sin 4 (/3/2) i cos , 



^co^'W^. (90) 



As £(/?) depends only on £+(/?) in the range (0, /?), it is 
possible to construct £(/?) in this range from an unbiased 
estimator of £+(/?) in the same range. By construction, 
£ (/3) — 5ft£_ (/3) will contain only B polarization in the mean. 

The window function for this method is simply -iKui 
(equation 57), so that 



<Cf > = 



(91) 



in the previous section, we can write -iKiy = 
-2Fi L -2W Ll i . Representative elements of the window 



As 

Ez 

functions -iKw are plotted in Fig. 6 for /3 max = 20° and 
Gaussian apodizing with half-width at half-maximum equal 
to /3max/2. They are well approximated by Gaussians with 
asymptotic normalisation given by the right-hand side of 
equation (82). For presentation purposes it is desirable to 
have window functions normalised to unity. We can enforce 
this by dividing the power spectrum reconstructed from 
equation (85) by y\, -2K H i . The exact normalisation is eas- 
ily computed from equation (65) by e.g. Gauss-Legendre in- 
tegration, and can be performed while inverting the correla- 
tion functions at negligible computational cost. The renor- 
malised window functions are also shown in Fig. 6, along 
with the mean recovered Cf ' , obtained from equation (85) 
with and without renormalisation by -2K H r . The renor- 
malised estimates agree very well in the mean with the true 
power spectra. Note also that since we have removed cross- 




fi (degrees) 

Figure 7. Azimuthally-symmetric weight function wp(n) 
adopted for the BICEP simulations. We chose to weight in pro- 
portion to the integration time per pixel. For white noise this is 
equivalent to weighting with the inverse of the noise variance. 



contamination, the recovered Cf are insensitive to large- 
scale power (from reionization) in E polarization. 



5.1 Application to BICEP 

As an application of our new estimator, we consider simu- 
lated maps for the BICEP experiment. BICEP is the first of 
a new generation of large bolometer arrays that are designed 
to target B-mode polarization. We used the experimental 
parameters taken from the BICEP homepage 6 . The survey 
will cover a polar-cap region of angular radius 18? 5, integrat- 
ing for a nominal 300 days. BICEP will be composed of 48 
polarization-sensitive bolometers (PSBs) operating at 100 
GHz with a resolution of 1° (full- width at half-maximum) 
and 48 at 150 GHz with 0?7 resolution. For our simula- 
tions we ignored the difference in beam size between the two 
channels taking the beam size to be 1°, so maps from each 
channel could be easily combined without introducing noise 



http://bicep.caltech.edu 
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Figure 6. Left: window functions _2-ft"ii' (top) and their rcnormaliscd counterparts (bottom) for /3 m ax = 20° and Gaussian apodization 
with half-width at half-maximum equal to /3 m ax/2. Right: mean recovered Cf with (crosses) and without (circles) renormalisation for 
the cosmological models with (bottom) and without (top) reionization. 
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Figure 8. The mean in 100 BICEP simulations of the flat band-powers for Cf (left; smoothed with a 1-degree beam) and Cf (right). 
Monte-Carlo error estimates from the 100 simulations are shown as boxes, centred on the average bandpowers from the simulations. The 
error bars are the theoretical approximation in equation (92). The solid lines are the input theoretical spectra, smoothed with the beam, 
for the cosmological model described in Sec. 1 with reionization at z = 6. Note they are not convolved with the band-power window 
function. 



correlations. We took each PSB to have an instantaneous 
sensitivity of 300 piK^/sec. 

We simulated 100 noisy CMB maps using a realistic 
map of the integration time per pixel based on the BICEP 
scanning strategy (see Fig. 7). The cosmology was that de- 
scribed in Sec. 1 (with reionization at z — 6). The pseudo- 
C\s were extracted with HEALPix at resolution iV s id = 512 
using the weight function shown in Fig. 7. This corresponds 
to inverse noise variance weighting for white noise in the 
time domain, and the azimuthal symmetry reflects the sym- 
metry of the proposed BICEP scan strategy. We generated 
a further suite of Monte-Carlo noise realisations which were 
used to remove the noise bias. The integral in equation (90) 
was performed with a cumulative Simpson rule, giving esti- 
mates of £ (/3) at the roots of a Legendre polynomial scaled 
to the angular range (0,31°). In principle we can estimate 
the correlation functions in the range (0, 37°), but very few 
pixel pairs contribute to the largest separation angles so the 
correlation functions are very noisy there. During the inte- 
gration to form £ (P) we constructed £+ (P) directly from the 



pseudo-C;s, with the noise bias removed, at points linearly 
spaced between the Legendre roots. (We used a nine-point 
Simpson rule.) 

We recovered the angular power spectra by evaluat- 
ing equation (85) with Gauss-Legendre quadrature. We 
adopted a Gaussian apodizing function with half-width at 
half- maximum equal to 18? 5. We further compressed our es- 
timates into flat band-powers with a Al = 35 thus removing 
much of the sensitivity to the choice of apodization. We ver- 
ified with Monte-Carlo simulations that Al = 35 is sufficient 
to remove any significant correlations between adjacent band 
powers. The mean band-powers for Cf and Cf (smoothed 
with the 1° beam) from 100 simulations are plotted in Fig. 8, 
along with ±<r error boxes estimated from the simulations. 
From the simulations we have verified that the method is 
unbiased (to within the standard error of the Monte- Carlo 
averages). We also compared the errors estimated from the 
100 simulations with the rule of thumb in equation (55). We 
found that to get good agreement with the Monte-Carlo er- 
rors on Cf it was necessary to refine equation (55) to take 
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account of the noise inhomogeneity and the compression to 
band powers more properly. We used the following approx- 
imation to the covariance of the recovered beam-smoothed 
Cfs (see Efstathiou 2003 for a derivation of this formula for 
the temperature anisotropies) : 

L M 



+ lyjCfC* ^2(wU 2 p ) L m(w 2 p )1 

M 




Here, ap(n) is the polarization noise variance per solid angle 
(sec equation 47), wp(h) is the polarization weight function, 
and e.g. (wpap)LM are the (spin-0) multipoles of the prod- 
uct w P (n)a'p(h). Equation (92) makes a number of approxi- 
mations: (i) it takes no account of the need to separate E and 
B, which is reasonable for E given that it dominates B, but 
will not be valid if extended to B when Cf jCf > (IA W ) 2 , 
where A w is the characteristic width of wp,i; (ii) it ignores 
the spin-2 nature of the polarization, which is acceptable at 
high and (iii) it only treats the inversion from pseudo- 
CiS to Cis approximately (i.e. divide by i02/ s k y ). We defer a 
full discussion of analytic approximations to the covariance 
of polarization power spectra to a future paper (Challinor 
& Chon in preparation), where we show how to generalise 
equation (92) to take account of E-B mixing on an apodized 
sky. Here, we simply note that the above assumptions should 
hold well for Cf in this application to BICEP. For the level 
of Cf in our assumed cosmology (r = 0.31), and given the 
smooth apodization of the edges of the survey region by 
wp(n), the application of equation (92) to B should serve 
as a useful first approximation to the errors. Note also that, 
for uniform noise, we recover equation (55) if we average into 
bands that are wide compared to the power spectrum of the 
square of the weight function, ^2 m \ {w 2 P )i m \ 2 / (21 + 1). The 
theoretical approximation to the errors in Fig. 8 are obtained 
by summing equation (92) over I and I' in a given band. 
These theoretical predictions agree well with the Monte- 
Carlo errors for E, and are in broad agreement for B. As 
expected, in the latter case the details of the E-B separation 
process that are ignored in our rough theoretical predictions 
are more critical. 

Our new estimator removes the cross-contamination be- 
tween E and B in the mean, however, the variance of an 
estimate of Cf or Cf contains a contribution from both E 
and B modes. To assess more carefully the level of cross- 
contamination due to the geometric effect of E-B mixing, 
we compute exactly the covariance of our decoupled power 
spectrum estimates in the absence of noise. Since the cross- 
contamination will be more significant for B than E, we 
concentrate on the former. We compute the error covari- 
ances first with E and B power retained, and then with 
only the B power. The latter calculation approximates the 
errors we would obtain if we separated the E and B modes 
at the level of the map prior to power spectrum estimation. 
The mechanics of the calculation are as follows. First, we 
compute the covariance of the polarization pseudo-C;s us- 
ing the techniques described by Hansen & Gorski (2003). 
This is only tractable beacuse of the azimuthal symmetry 
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Figure 9. The variance of an estimate of Cf obtained by equa- 
tion (93) in the absence of noise with the BICEP sky coverage. 
The boxes are the full error, including the cross-contribution from 
E modes. The bars are the variance obtained by setting Cf to 
zero, and thus are representative of the errors that would be 
achieved if B were separated from E at the level of the map. 



of the BICEP scanning strategy. We then linearly transform 
the pseudo-C; covariance to that of the decoupled estimates 
using the fact that our power spectrum estimation method 
is linear in the pseudo-C;s, i.e. 

cov{C?,6J,)= M? L x 'cov(C?,CZ;)MZZ;, (93) 

L,L',X',Y' 

where X, X' , Y and Y' run over E and B, and the coupling 
matrices relate the decoupled power spectrum estimates to 
the pseudo-C;s: 

I'X' 

The columns of the coupling matrices M x , x are conve- 
niently extracted by setting all of the pseudo-C; s to zero ex- 
cept for C x in our power spectrum code, and then reading 
out the recovered C x . The matrix is symmetric on the in- 
dices X and X' . Figure 9 summarises our results obtained for 
the B-mode power spectrum. As in Fig. 8 we compress our 
results into flat band-powers. We see that, in the noise-free 
case, the errors we obtain using the BICEP weight function 
wp(n) are not dominated by the cross-contribution from 
i?-mode polarization; this accounts for approximately only 
20-per cent of the error budget. Of course, a lower level of B- 
mode polarization would increase the relative significance of 
the cross-contribution. Ultimately, this cross-contamination 
may make our method unsuitable for future, high-sensitivity 
B-mode experiments surveying small regions if the tensor 
amplitude is too low. We leave quantification of this state- 
ment to a future paper (Challinor & Chon in preparation), 
where we show how to estimate the cross-contribution to the 
variance for non-symmetric weight functions. 

Since full covariance information is available from equa- 
tion (93), the level of correlation between adjacent band- 
powers can be calculated directly. For the BICEP scan strat- 
egy, M X , Y oscillates with a full- width at half-maximum of 
approximately 12. Hence, for the chosen bin width Al — 35, 
the correlation between adjacent bins is negligible. 



© 2003 RAS, MNRAS 000, 1-14 



Fast polarization power spectrum estimation 13 



6 CONCLUSION 

We presented a fast and unbiased method to extract CMB 
polarization power spectra from large maps via the two- 
point correlation functions. The method, which generalises 
that of Szapudi et al. (2001) to polarization, can be sum- 
marised as follows. First, we compute unbiased estimates of 
the three (complex) polarization auto- and cross-correlation 
functions at the roots of a Legendre polynomial from pseudo- 
Cis of heuristically- weighted maps. The estimates of the cor- 
relation functions can be computed in O(N^) operations 
using fast spherical transforms. If the correlation functions 
can be estimated for all angular separations, the power spec- 
tra can be accurately recovered with Gauss-Legendre inte- 
gration. In this case, the method is unbiased: the theoreti- 
cal window function is a Kronecker delta. Further compres- 
sion to band-powers can then be made, and the resulting 
theoretical window functions would be top-hat functions. If 
the correlation functions cannot be estimated for all angular 
separations, due to limitations of sky coverage, we showed 
that significant E-B mixing can occur. In particular, large- 
scale E power (due to reionization) can be aliased into B 
on all scales. Although E-B mixing does not present a fun- 
damental problem for parameter extraction, it does compli- 
cate the interpretation (and presentation) of the recovered 
power spectra. For this reason, we proposed a new estima- 
tor, extending earlier work by Crittenden et al. (2002), that 
removes E-B mixing exactly in the mean when working with 
incomplete correlation functions. Note that this is not the 
case for regularised inversions of the pseudo-C;s in harmonic 
space (e.g. by working with pseudo band powers as in a 
polarized extension of MASTER; Hivon et al. 2002). The 
new estimator requires one further numerical integration of 
the estimated correlation function £+{tp) to obtain functions 
that contain only E or B power in the mean. Using e.g. a 
cumulative Simpson rule, these functions can be estimated 
accurately at the roots of a Legendre polynomial and in- 
verted to power spectra with Gauss-Legendre quadrature. 
The increase in computational effort is minimal, and the 
theoretical window functions that result do not couple E 
and B power in the mean by construction. Fourier ringing 
in the estimates can be safely controlled by apodizing the 
integral transforms (or by compressing into band-powers), 
without introducing any E-B mixing. 

An essential part of our method (and indeed any 
quadratic method) is being able to remove the mean noise 
contribution from the correlation functions (i.e. to remove 
the bias due to the noise). For general noise properties we 
must resort to Monte-Carlo evaluation of the mean over an 
ensemble of pure noise realisations. The method presented 
here is thus dependent on being able to simulate noise maps 
efficiently. Error estimation on the recovered C;S must also 
generally proceed by Monte-Carlo evaluation. The O(N^) 
scaling of our method makes this a realistic proposition, even 
for mega-pixel maps. 

We applied our methods to simulations of a large-area 
survey, with parameters similar to Planck, and also to the 
BICEP experiment which will cover 3 per cent of the sky. 
In both cases we obtained errors in line with theoretical 
expectations. Although our algorithm in its present form 
is already practical for analysing mega-pixel CMB maps, 
further work is required to assess the optimality of the un- 



derlying methods. In particular, comparison with current 
(brute-force) maximum-likelihood codes should be possible 
for low-resolution simulations; comparison at higher reso- 
lution must await further algorithmic development of the 
likelihood codes. Another issue worth investigating further 
is the impact of the choice of pixel weighting on the cosmic 
variance contribution from e.g. Cf to the recovered Cf . Al- 
though we have separated E and B in the mean, this does 
not guarantee that in a single realisation there is no recov- 
ered B power due to E modes in the signal. We have shown 
that this geometric effect does not dominate the error bud- 
get for Cf for the BICEP weight function, in the absence of 
noise, for a tensor-to-scalar ratio of r — 0.31. However, ulti- 
mately this cross-contribution may limit the applicability of 
our method if the tensor amplitude is low enough (Challi- 
nor & Chon in preparation). To eliminate this undesirable 
contribution to the variance would require separating E and 
B at the level of the map, e.g. Lewis et al. (2002), prior to 
performing power spectrum estimation. 
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APPENDIX A: A(ip) FOR 
UNIFORM ALLY- WEIGHTED, 
AZIMUTHALLY-SYMMETRIC PATCHES 

In the special case that the analysis is performed over an 
azimuthally-symmetric part of the sky with uniform pixel 
weighting w(n) = 1, the correlation function normalisation 
A(ip) can be evaluated analytically. (We can drop the sub- 
scripts P, T and X in this case since there is no distinction 
between the normalisations.) 

We begin by considering a polar-cap region with an- 
gular radius a, and assume that a < 7r/2. The integral in 
equation (31) then evaluates to give 



4tt 



AW 



cos 1 



2sin 2 (^/2) 



— cos a cos 1 



2tan 2 0/>/2) 
tan 2 a 



(Al) 



for ip 5j 2a, and zero otherwise. Note that 1/A(Q) = 47r 2 (l — 
cos a) = 87r 2 / sky as required by equation (53). Note also 
that 1/A(ip) goes continuously to zero as ip — > 2a. We can 
now construct the other important case of a Galactic cut 
by symmetry. For a symmetric cut subtending an angle a c , 
provided that a c < 7r/2, there are pixel pairs at all angular 
separations and 1/A(ip) is non-zero for all ip € (0, it). If we 
denote the normalisation for a polar-cap region of angular 
radius (it — a c )/2 by A*(tp), then we find 

+ A,(l-ti) ac<^<vr-Q c , (A2) 
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A„(y>) A,(-n-i>) 
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A* (7T — 



7T — Ct c < Ip < 7T. 
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